
options(stringsAsFactors=F)
library(data.table)
library(parallel)
library(Seurat)

### testis01 ###
all_matrix_01 <- Read10X(data.dir = "/public/user/yl2019/20211229_testis_singleCell_V2/20220107_exp_atac_res_V1/outs/filtered_feature_bc_matrix")
exp_matrix_01 <- all_matrix_01[["Gene Expression"]]  ## 36601 gene, 20000 cell ##
rna_01 <- CreateSeuratObject(counts = exp_matrix_01, assay = 'RNA', min.cells = 0, min.features = 0, project = '10x_RNA')  

### testis02 ###
all_matrix_02 <- Read10X(data.dir = "/public/user/yl2019/20220926_testis/testis_02/filtered_feature_bc_matrix")
exp_matrix_02 <- all_matrix_02[["Gene Expression"]]  ## 36601 gene, 9815 cell ##
rna_02 <- CreateSeuratObject(counts = exp_matrix_02, assay = 'RNA', min.cells = 0, min.features = 0, project = '10x_RNA')  

### testis03 ###
all_matrix_03 <- Read10X(data.dir = "/public/user/yl2019/20220926_testis/testis_03/filtered_feature_bc_matrix")
exp_matrix_03 <- all_matrix_03[["Gene Expression"]]  ## 36601 gene, 11725 cell ##
rna_03 <- CreateSeuratObject(counts = exp_matrix_03, assay = 'RNA', min.cells = 0, min.features = 0, project = '10x_RNA')  


### 数据合并 ###
rna_all <- merge(rna_01, y = c(rna_02, rna_03), add.cell.ids = c("testis01", "testis02", "testis03"), project = "testis")  ## 41540 cell ##
mito.features <- grep(pattern = "^MT-", x = rownames(x = rna_all), value = TRUE)
percent.mito <- Matrix::colSums(x = GetAssayData(object = rna_all, slot = 'counts')[mito.features, ]) / Matrix::colSums(x = GetAssayData(object = rna_all, slot = 'counts'))
rna_all$percent.mito <- percent.mito
pdf(file='/public/user/yl2019/20221201_testis_singlem_upd/ATAC_res/correlation_20221214/quality.pdf',width=8,height=6)
VlnPlot(object = rna_all, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), group.by="batch",pt.size = 0)
dev.off()

rna_all@meta.data$cell <- rownames(rna_all@meta.data)
rna_all@meta.data$batch <- sub("(\\w+)_(\\w+)-(\\w+)","\\1",rna_all@meta.data$cell)

## 质控标准画图 ##

pdf(file='/public/user/yl2019/20221201_testis_singlem_upd/RNA_res/quality.pdf',width=8,height=6)
VlnPlot(object = rna_all, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
dev.off()


rna_all_QC3 <- subset(x = rna_all, subset = nFeature_RNA > 500 & nFeature_RNA < 10000 & nCount_RNA > 500 & nCount_RNA < 50000 & percent.mito < 0.25)  
rna_all_QC3 <- rna_all_QC3[,cell_QC3passf$cell]

rna_all_QC3_magic <- NormalizeData(object = rna_all_QC3)
rna_all_QC3_magic <- ScaleData(object = rna_all_QC3_magic)
rna_all_QC3_magic <- magic(rna_all_QC3_magic, knn=10, t=6, npca=20)

# rna_all_QC3_magic@active.assay = "RNA"
testis.list <- SplitObject(rna_all_QC3_magic, split.by = "batch")

# normalize and identify variable features for each dataset independently
testis.list <- lapply(X = testis.list, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})

# select features that are repeatedly variable across datasets for integration
testis_features <- SelectIntegrationFeatures(object.list = testis.list)
testis_anchors <- FindIntegrationAnchors(object.list = testis.list, anchor.features = testis_features)
testis_combined <- IntegrateData(anchorset = testis_anchors)

DefaultAssay(testis_combined) <- "integrated"
# Run the standard workflow for visualization and clustering
testis_combined <- ScaleData(testis_combined, verbose = FALSE)
testis_combined <- RunPCA(testis_combined, npcs = 100, verbose = FALSE)
testis_combined <- RunUMAP(testis_combined, reduction = "pca", dims = 1:30)

# testis_combined <- RunHarmony(testis_combined, "batch" , plot_convergence = FALSE, dims.use = 1:30, reduction.save = "harmony")
# testis_combined <- RunUMAP(testis_combined, reduction = "harmony", dims = 1:30)

testis_combined <- FindNeighbors(testis_combined, reduction = "pca", dims = 1:30)
# testis_combined <- FindNeighbors(testis_combined, reduction = "harmony",dims = 1:30)
testis_combined <- FindClusters(testis_combined, resolution = 0.3)
testis_combined <- RunUMAP(object=testis_combined, dims = 1:30)

## 画聚类图 ##
DefaultAssay(testis_combined) <- "MAGIC_RNA"
pdf(file='/public/user/yl2019/20221201_testis_singlem_upd/RNA_res/QC3_magic3_res/QC3_magic3_cluster.pdf',width=8,heigh=6)
DimPlot(testis_combined, reduction = "umap", label = TRUE)
dev.off()

## marker基因 ##
pdf(file='/public/user/yl2019/20221201_testis_singlem_upd/RNA_res/QC3_magic3_res/QC3_magic3_cluster_marker.pdf',width=8,heigh=6)
FeaturePlot(testis_combined, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP","CD8A"))  ## 这个基因名不是睾丸的marker基因,你要画的话需要改一下
dev.off()
